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ABSTRACT 

We generate simulations of the CMB temperature field as observed by the WMAP satellite, taking 
into account the detailed shape of the asymmetric beams and scanning strategy of the experiment, 
and use these to re-estimate the WMAP beam transfer functions. This method avoids the need of 
artificially symmetrizing the beams, as done in the baseline WMAP approach, and instead measures 
the total convolution effect by direct simulation. We find noticeable differences with respect to the 
nominal transfer functions. For instance, the nominal VI beam under-estimates the full beam convo- 
lution by ~ 0.5% a.t £ — 500 and ~ 1.0% at £ = 800. Similar differences are seen for other DA's. This 
in turn implies that the high-£ power spectrum is biased low by 1 — 2%, effectively tilting the spectrum 
slightly. Re-estimating cosmological parameters we find that the spectral index of scalar perturbations 
is Us = 0.969 ±0.014 after correcting for this effect, corresponding to a positive shift of O.Scr compared 
to the previously released WMAP results. Our CMB sky simulations are made publicly available, and 
can be used for general studies of asymmetric beam effects in the WMAP data. 
Subject headings: cosmic microwave background — cosmology: observations — methods: statistical 



1. INTRODUCTION 

Without doubt, the angular CMB power spectrum 
is today our single most important source of cosmo- 
logical information. Perhaps the most striking demon- 
stration of this fact t o date is the WMAP expe riment , 
(jBennett et al.l 120031 : iHinshaw etall l2007l . |2009[) which 
has allowed cosmologists to put unprecedented con- 
straints on all main cosmological parameters, as well as 
ruling out vast regions of the possible model spaces. Sim- 
ilarly, in only a few years from now Planck will finally 
provide the definitive measurements of the temperature 
power spectrum, as well as polarization spectra with un- 
precedented accuracy. This will certainly lead to similar 
advances in our knowledge about the history of our uni- 
verse. 

Each of these experiments observes the CMB field by 
scanning the sky with an instrumental beam of finite res- 
olution. This operation effectively corresponds to averag- 
ing over beam-sized angular scales, and is expressed tech- 
nically either in pixel space by a convolution of the beam 
with the underlying sky, or in harmonic space by a mul- 
tiplication of the two corresponding sets of harmonic ex- 
pansion coefficients. For simplicity, the harmonic space 
expansion of the beam is typically expressed in terms of 
Legendre coefficients of an (azimuthally symmetric) ef- 
fective beam response. This function is often called "the 
beam transfer function", bi. 

Before it is possible to make unbiased cosmological in- 
ferences based on the CMB power spectrum, it is of crit- 



s: |i.k.wehus@fys.uio.no| 

3: lottyQcaltech.edul 



Electronic address: 

Electronic address: lottyCScaltech 

Electronic address: h.k.k.eriksen@astr o.ulo.nol 
Electronic address: leuat@irio.co.uk 

^ Department of Physics, University of Oslo, P.O. Box 1048 Blin- 
dern, N-0316 Oslo, Norway 

^ California Institute of Technology, Pasadena, CA 91125 
^ Institute of Theoretical Astrophysics, University of Oslo, P.O. 
Box 1029 Bhndern, N-0315 Oslo, Norway 

Centre of Mathematics for Applications, University of Oslo, 
P.O. Box 1053 Blindern, N-0316 Oslo, Norway 



ical importance to know the beam transfer function to 
high precision, as an error in the beam function translates 
into a direct bias in the estimated power spectrum. This 
in turn requires detailed knowledge about the beam re- 
sponse function on the sky for each experiment. For a full 
description of t he WMAP beam e s timation process and 
fina l model, see iPage et al.l (|2003f l. iJarosik et all (|2007l ) 
and lHill et all f^Om . 

The impact of asymmetric beams may also be impor- 
tant for applications other than power spectrum esti- 
mation. One example of special interest to us is the 
assessment of non-Gaussia nity and violatio n of statisti- 
cal isotropy. Specifically, lAckerman et al.l |2007,) con- 
sidered a model based on violation of rotational invari- 
ance in the early universe, and derived explicit paramet- 
ric expressions for the correspon ding observational signa- 
ture. Then, in a follow-up paper iGroeneboom &: EriksenI 
()2009() analysed the 5-year WMAP data with respect to 
this model and, most surprisingly, found a detection at 
the 3.8(7 confidence level. Given that this was a most un- 
expected result, several questions concerning systematic 
errors in the WMAP data were considered, in particular 
those due to residual foregrounds, correlated noise and 
asymmetric beams. However, it was shown in the same 
paper that neither foregrounds nor correlated noise were 
viable explanations, while the question of asymmetric 
beams was left unanswered, due to lack of proper sim- 
ulation machinery. This question provided our initial 
motivation for considering the problems studied in this 
paper. 

The starting point for tackling the asymmetric beam 
problem for WMAP is a set of beam maps released by the 
WMAP team, two for each differencing assembly (DA), 
denoted A and B, respectively. These maps were derived 
by observing Jupiter for extended periods of time. Then, 
in order to derive the proper beam transfer functions, the 
WMAP team adopted a computationally fast and con- 
venient approach: They first symmetrized the effective 
beam for each DA, collapsing the information in the A 
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and B sides into one common function, and then com- 
puted the Legendre transform of the corresponding radial 
profile. However, for this to be an accurate approxima- 
tion, one must on the one hand assume that the beams 
on the two sides are very similar, and on the other hand 
either assume that both beams are intrinsically circularly 
symmetric, or that all pixels on the sky are observed from 
all angles an equal number of times due to the scanning 
strategy. In reality none of these conditions are met, and 
one may therefore ask whether there might be any resid- 
ual effect due to the combination of an asymmetric beam 
and anisotropic scanning in the WMAP beam functions. 

This probl ern was addressed analytically by 

iHinshaw et al.l (|2007f ). who derived an approximate 
expression for the expected power spectrum bias due to 
asymmetric beams in the WMAP data. Their conclusion 
was that such effects were <1% for the 3-year WMAP 
data. 

In this paper, we revisit the question of asymmet- 
ric beams in WMAP with two main goals. First, we 
seek to estimate the effective beam transfer functions for 
each WMAP DA, taking into account the full details 
of the asymmetric beams and specifics of the WMAP 
scanning strategy by direct simulation. This way, we 
check whether the analytic approximations presented by 
[Hinshaw ct al. ( 200'i^ are valid. Second, we want to pro- 
duce a set of high-fidelity simulated CMB sky maps, with 
beam properties as close as possible to those observed 
by WMAP, that can later be used for general studies of 
asymmetric beam effects in WMAP. 

2. PIPELINE OVERVIEW 

In this section we summarize the methods and algo- 
rithms used in this paper. Note that none of the indi- 
vidual steps described below are original to this paper, 
and only the main ideas will therefore be discussed in the 
following. 

We begin by defining our notation. We will be estimat- 
ing the product of the WMAP beam transfer function, bi, 
and pixel window, pe, by direct simulation. This product 
is denoted (3i — bipi. Given this function, the combined 
effect on a sky map, T{n), of convolution by an instru- 
mental beam and averaging over finite-sized pixels may 
be approximated in harmonic space as 

£=0 m=~l 

where Yfm(ri) are the usual spherical harmonics. 
The angular power spectrum of T is given by 

1 ^ 

m=-£ 

while the power spectrum of the true underlying CMB 
map, sin), is 

1 ^ 

The effect of the beam convolution and pixel averaging 
on the power spectrum is therefore simply given by a 
multiplication with 



The overall approach for estimating [3^ used in this pa- 
per may be summarized by the following steps: First, we 
simulate time-ordered data (TOD) for each DA, taking 
into account both the detailed beam maps of WMAP 
and the exact orientation of the spacecraft at each point 
in time. We then produce a sky map from this TOD. 
Finally we compute the square root of the ratio between 
the output and the input power spectra, which becomes 
our estimate of (3g. 

Note that in this paper we are only concerned with the 
effect of asymmetric beams, not other systematic effects 
such as instrumental noise. All following discussions will 
therefore assume noiseless observations. 

2.1. Simulation of time-ordered data 

Our first step is to simulate a reference CMB sky re- 
alization, s, given an angular temperature power spec- 
trum, C^''°°'^. This can be achieved with a standard code 
such as "anafast", which is available in the HEALPix^ 
software package. Note that this map should not be 
smoothed with either an instrumental beam or a pixel 
window; adding these effects is the task of the following 
pipeline. Explicitly, the input reference map should sim- 
ply be pure spherical harmonic modes projected onto a 
set of pixel centers. 

Next, we need to be able to convolve this map with 
a given beam map at arbitrary positions and orienta- 
tions on the sphere. In this paper we do this by brute- 
force integration in pixel space. For an alternative fast 
Fourier space base d appr oach to the same problem, see 
iWandelt fc Cfekil (|200l . 

We define p to be a unit vector pointing towards the 
beam center, and specify its position on the sphere using 
longitude and co-latitude {(p, 9) . We further define -0 to 
be the angle between some fixed reference direction in 
the beam map and the local meridian. The value of the 
beam map at position h = ((/)', 0'), which in principle 
is non-zero over the full sky, is denoted b{(j)' ,9';4>,0,tp). 
With these definitions, the desired convolution may be 
written as 

T(</.,0,^)=/ s{cj,',e')b{<j>',9';cp,e,^p)dn'. (4) 

Ji-rr 

Computationally speaking, we approximate this inte- 
gral as a direct sum over HEALPix pixels, which all 
have equal area, with the product s ■ b being evaluated 
at HEALPix pixel centers. To make these calculations 
computationally feasible, we assume that the beam is 
zero beyond some distance from the beam center (rang- 
ing between 3.5 and 7° for the WMAP channels), and 
thus only include the main lobe in the following analysis. 
While the WMAP beam maps are provided as pixelized 
maps, we need to know the beam values at arbitrary po- 
sitions (ie., HEALPix pixel centers). We solve this by 
computing a 2D spline for each beam map, enabling us 
to interpolate to arbitrary positions. For a review on one 
specific method for fast (and local) 2D spline evaluations, 
see Appendix [Bl 

WMAP is a differential experiment, and measures at 
each point in time the difference between the signals re- 
ceived by two different detectors, denoted A and B. The 

® http://hcalpix.jpl.nasa.gov 
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full set of time-ordered WMAP data may therefore be 
written as 



d,{^)^T^{^)-T^{^), 



(5) 



where x = {Kl, Kal, Ql-2, Vl-2, Wl-4} is a DA label, 
and j is a time index, and for each detector a short- 
hand for {(f), 9, ip). This equation may be written in the 
following matrix form. 



AT, 



(6) 



where we have introduced an A^tod x A'pix pointing matrix 
A. This matrix contains two numbers per row; 1 in the 
column hit by the center of beam A at time i, and -1 in 
the column hit by the center of beam B. 

The remaining problem is to determine the position 
and orientation of each detector at each time step. This 
information has been made publicly available by the 
WMAP team on LAMBDA^, and consists of a large set 
of pointing files together with useful IDL routines for 
extracting the desired information. 

2.2. Map making with differential data 

For the ma p making step we adopt the algorithm de- 
veloped bv Wright et al.l (llQQGl). which was used in the 1 - 
and 3-year WMAP pipelines ("Hinsha w et al.ll2003l l2007f). 
Here we only summarize the essential algebra, and out- 
line the algorithm. 

Our goal is to establish an unbiased and, preferably, 
optimal estimate of the (smoothed) sky signal, T, given 
a set of differential TOD values, d. For noiseless data, 
the maximum likelihood estimator is simply 



(A*A)-iA*d. 



(7) 



For high-resolution sky maps, this equation involves 
an inverse of a large matrix and cannot be solved ex- 
plicitly. Instead, one often resorts to iterative methods 
such as Conjugate Gr adients, or, for differe ntial data, the 
method developed bv lWright et al.l (|1996D . 

We present the iterative differential map maker in a 
simple manner: Define D to be the diagonal matrix that 
counts the number of hits Aobs(p) per pixel p on the diag- 
onal, and Oj and hi to be the pixels hit by side A and side 
B at time i, respectively. Suppose that we already have 
established some estimate for the solution, T-'. (Note 
that this can be zero.) Then the iterative scheme 



(8) 



will converge to the true solution: If = T, then d = 
AT-' , and the second term on the right hand side is zero. 
This algorithm is implemented by the following scheme: 



rpj+l 

V 



(9) 



This algo r ithm was originally presented by 
iWright et all (|l99l . The only new feature intro- 
duced here is the cho i ce of s tarting point. In the original 
paper, IWright et al.l (|1996D initialized the iterations at 
the DMR dipole, since their test simulation included 
a CMB dipole term. However, for a given scanning 
strategy, there will often be some large-scale modes that 

® http://lambda.gsfc.nasa.gov 
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Fig. 1. — Top panel: Comparison between the transfer functions, 
/Bi, for a Gaussian beam of 20' FWHM. This was computed from 
a pixelized beam map and with the WMAP VI scanning strategy 
(red line), and alternatively, by using the well-known analytic ex- 
pression for the Legendre transform of a Gaussian beam (Equation 
I12| l and isotropized HEALPix pixel window (black dashed line). 
The vertical dotted line indicates the multipole moment, ^hybrid, 
at which /3f = 0.15. Bottom panel: The ratio between the trans- 
fer functions in the top panel. Note the excellent agreement up to 
£ ~ 800, after which the differences in pixel window approximations 
becomes visible. 
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Fig. 2. — Comparison of the VI WMAP beam transfer functions 
obtained from six (dashed line) and twelve (solid line) months of 
observations. Note that the behaviour at low and intermediate ^'s 
is similar, whereas the pixel window-induced increase is seen earlier 
in £ for the six month data set than for the twelve month data set. 

are less well sampled than others. For instance, for 
the WMAP s trategy £ = 5 is mo re problematic than 
other modes (jHinshaw et al.|[2003f ). This leads to slow 
convergence with the above scheme for this mode. 

We therefore choose a different approach: Before solv- 
ing for the high-resolution map by iterations, we solve 
Equation[7|by brute-force at low resolution. For the cases 
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TABLE 1 
Summary of DA parameters 



DA 


FWHM'» 


Radius'" 


A^sidc 




■^hybrid 


^ samples 






(arcmin) 


(degrees) 








(10*^) 


(mK) 


Kl 


53 


7.0 


512 


750 


318 


2.5 




Kal 


40 


5.5 


512 


850 


411 


2.5 




Ql 


31 


5.0 


512 


1100 


522 


3.1 


78.2 


Q2 


31 


5.0 


512 


1100 


515 


3.1 


74.2 


VI 


21 


4.0 


1024 


1500 


789 


4.1 


99.0 


V2 


21 


4.0 


1024 


1500 


779 


4.1 


88.2 


Wl 


13 


3.5 


1024 


1700 


1164 


6.2 


143.8 


W2 


13 


3.5 


1024 


1700 


1148 


6.2 


159.7 


W3 


13 


3.5 


1024 


1700 


1162 


6.2 


168.5 


W4 


13 


3.5 


1024 


1700 


1169 


6.2 


164.4 



^ Effective symmetrized beam size.'' Radius used for pixelized beam con- 
volutions. See Hill ct al] I I2009I) for details. Average full-sky RMS values 
evaluated at A^sidc = 512. 



considered later in this paper, we choose a HEALPix 
resolution of A^sidc = 16 for this purpose. With 3072 
pixels, about 30 seconds are needed to solve this sys- 
tem by singular value decomposition. (Note that the 
monopole is arbitrary for differential measurements, and 
one must therefore use an eigenvalue decomposition type 
algorithm to solve the system.) The improvement in con- 
vergence speed due to this choice of initial guess is ex- 
plicitly demonstrated in Appendix [Xl 

Our convergence criterion is chosen such that the RMS 
difference between two consecutive iterations must be 
less than 0.05 fxK. We have verified that this leads to 
errors of less than 0.1 fiK in the final solution, of which 
far most is due to a residual dipolc. This is typically 
achieved with 30-50 iterations, although some converge 
already after 20-30 iterations and a few after 70 or more 
iterations. 

At first glance, the fact that the final residuals are as 
small as 0.1 fiK for an RMS stopping criterion as large 
as 0.05 /xK may seem surprising. However, this is ex- 
plained by the fact that the iterative solution obtained 
by Equation [8] often alternates between high and low val- 
ues about the true answer. This suggests that a further 
improvement to the algorithm may be possible: Faster 
convergence may perhaps be obtained by computing the 
average oftwo consecutive iterations, T = {T^ +T^+'^)/2, 
as the map estimate for iteration j + 2. However, the 
computational resources spent during map making is by 
far sub-dominant compared to the TOD simulation, and 
we have therefore not yet implemented this step in our 
codes. 

2.3. Estimation of hybrid beam transfer functions 

As described in the introduction to this section, we 
estimate the transfer function by the square root of the 
ratio between the power spectra of the convolved map 
and the input map 



/3| 



(10) 



However, as noted above, this function describes both 
the effect from instrumental beam smoothing and aver- 
aging over pixels. In the present paper we are concerned 
mostly with the former of these, which has a stronger im- 
pact on large to intermediate scales. This is because the 
beam component is largely independent of total obser- 
vation time, assuming at least one year of observations 



for WMAP, whereas the pixel averaging component de- 
pends strongly on total observation time, or the average 
number of samples per pixel. The latter therefore evolves 
much more strongly with time than the former, as will 
be explicitly demonstrated later. 

In order to provide transfer functions that are valid for 
long observation periods (e.g., 5 or 7 years of WMAP 
observations), we choose to construct a hybrid transfer 
function, 




for £< i 



hybrid 



for e > ^hybrid 



(11) 



Here ^wmap ^j^^ nominal symmetrized transfer func- 
tion published by the WMAP team, is the (uniformly 
averaged) HEALPix pixel window, and ^hybrid is some 
transition multipole. In other words, we adopt our own 
direct estimate of the transfer function up to i^hybrid, but 
the symmetrized, asymptotically uniform and properly 
scaled WMAP transfer function at higher multipoles. 

Note that this issue is of minor importance in terms of 
cosmological interpretation, i.e., angular power spectrum 
and cosmological parameters, because the transition typ- 
ically takes place in the noise dominated high-^ regime. 
The effect of the anisotropic pixel window is therefore 
largely suppressed. In the present paper, we therefore 
choose to focus on the beam dominated region, and leave 
a detailed study of the pixel window to a future paper. 
See Section |4] for a detailed discussion of this issue. 

Finally, because we only generate a relatively small 
number of simulations in this paper, there is considerable 
Monte Carlo scatter in our estimated transfer functions 
on an i-hy-£ basis. To reduce this Monte Carlo noise, 
we smooth all transfer functi ons using the smooth splin e 
formalism described by, e.g.. lGreen fc SilvermanI (|1994l ). 

3. DATA AND SIMULATIONS 

All data products used in this study are provided by 
the WMAP team on LAMBDA as part of their 5-year 
data release. However, the calculations performed here 
are computationally extremely demanding, and we there- 
fore include only roughly one year worth of data in our 
calculations. To be precise, we include the period be- 
tween July 10th 2001 and August 2nd 2002, except for 
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Fig. 3. — Comparison between transfer functions derived in this paper (red lines) to the nominal WMAP transfer functions (black dashed 
lines). The transition multipole, ^hybrid is marked by dotted vertical lines. 
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three days with missing data, for a total of 383 days'". 

We consider all 10 WMAP DAs in our calculations, 
which are denoted, in order from low to high frequen- 
cies, Kl (23 GHz), Kal (33 GHz), Ql-2 (41 GHz), Vl-2 
(61 GHz) and Wl-4 (94 GHz), respectively Their reso- 
lutions range from 53' FWHM at K-band to 13' FWHM 
at W-band. Because of this large range in resolution, we 
specify the pixel resolution and harmonic space range for 
each case separately. For instance, K-band is pixelized at 
-^side = 512, and includes multipoles up to ^,nax = 750 
(the highest multipolc present in the transfer function 
provided by the WMAP team), while the W-band is pix- 
elized at A'^sido — 1024, and includes multipoles up to 
^max = 1700. A full summary of all relevant parameters 
for each DA is given in Table [T] 

Note that the listed noise RMS values are only used for 
estimating the power spectrum weights in Section [S] For 
simplicity we have adopted the official RMS values for 
the foreground-reduced 5-year WMAP maps here, but 
note that there is a ^ 1% bias in some of these values 
(|Groeneboom et aI]|2009D . However, this has no signifi- 
cant impact on the results presented in this paper. 

The beam maps for each DA are provided in the form 
of pixelized maps, and separately for side A and B. 
Each beam map contains non-zero values inside a radius 
around the beam center which is specified for each DA. 
For instance, the K-band radius is 7°, and the W-band 
radius is 3? 5. When evaluating the convolution defined 
in Equation m we include all pixels inside this radius. 

The pixel size of the beam maps is 2.4', which over- 
samples even the W-band beams. Based on these high- 
resolution maps we precompute all coefficients of the cor- 
responding bi-cubic spline (see Appendix [B|) , which al- 
lows us to very quickly interpolate at arbitrary positions 
in the beam map with high accuracy. 

Each beam is normalized by convolving a map constant 
equal to 1 at 1000 random random positions and orien- 
tations, and demanding that the average of the resulting 
1000 values equals unity. With the 2D spline interpola- 
tion scheme described in Appendix |B] the random uncer- 
tainties on the normalization due to beam position and 
orientation are ^ 0.2%. For comparison, directly reading 
off pixel values from the beam maps without interpola- 
tion leads to variations in the normalization at the ^ 2% 
level. 

For our base CMB reference sky set, we draw ten ran- 
dom Gaussian realizations from the the best-fit ACDM 
power spectrum derived fro m the 5-year WMAP data 
alone (|Komatsu et al.ll2009f ). These maps are generated 
at both (iVsidG = 512 and iVsidc — 1024 using the same 
seeds, and include neither an instrumental beam nor a 
pixel window; they are simply spherical harmonic modes 
projected onto the HEALPix pixel centers. All ten real- 
izations are processed for all ten DAs, such that the re- 
sulting simulations may be used for multifrequency anal- 
ysis, if so desired. 

As noted above, the computational requirements for 
the analyses presented here are demanding. The CPU 

Our original intention was to include precisely one year of ob- 
servations in our analysis, and therefore we processed 365 WMAP 
pointing files. However, we noticed after the calculations were com- 
pleted that some of the pointing files contained slightly more than 
one day's worth of data, such that a total of 383 days was in fact 
included. 



time for processing a single W-band DA is ~ 4000 hours, 
and the total disk usage for the entire project is '^ITB. 
For comparison, the corresponding map making step re- 
quires ~ 60 CPU hours, and is thus completely sub- 
dominant the TOD simulation. 

4. COMPARISON WITH ANALYTIC CASE 

In order to test our pipeline and understand its out- 
puts, we start by considering a perfect Gaussian beam. 
This case is treated in two different ways: First, we con- 
volve a CMB realization directly in harmonic space (as 
defined by Equation [T]) with a crfwhm = 20' FWHM ana- 
lytic Gaussian beam and the appropriate HEALPix pixel 
window, pi for A^sidc = 1024. The combined transfer 
function for this case reads 

f3f = e-^^(^+i)'^Vi, (12) 

where a = crf„hm/V8 ln2, and crtwhm is expressed in ra- 
dians. 

Second, we map out a corresponding two-dimensional 
Gaussian in pixel space over a grid of 2.4' pixels, the 
same resolution as the WMAP beam maps. We then 
input this into our simulation pipeline together with the 
same CMB realization used for the analytic convolution, 
and with the VI channel pointing sequence. From the 
resulting brute-force convolved map, we then obtain the 
effective transfer function, (3e, as described in Section [2T3l 

This function is plotted in the top panel of Figure [1] 
together with the product of the analytic Gaussian beam 
and the HEALPix pixel window. The ratio of the two 
effective functions is shown in the lower panel. 

From this figure it is clear that the agreement between 
the two approaches is excellent up to « 800. At higher 
^'s, however, the ratio increases rapidly, indicating that 
the analytic approach smooths more than the brute- force 
approach. This is due to the different definitions of the 
pixel windows in the two cases: In the HEALPix case, 
the pixel window is defined as an effective average both 
over each pixel and over the full sky. In other words, it 
assumes that all points have been observed an equal (and 
infinite) number of times. 

However, this is not the case for a real experiment 
which scans the sky for a finite length of time. As a con- 
sequence, each pixel is observed only a relatively small 
number of times, and this leads effectively to less smooth- 
ing. In the extreme case of only one observation per pixel, 
there would be no pixel averaging at all. Now, the av- 
erage pixel window would formally equal unity. On the 
other hand, the random realization-specific uncertainties 
would be very large, and the pixel window as such would 
have zero predictive power. 

In practice, one is well advised not to consider scales 
smaller than those that are properly oversampled by the 
scanning strategy. In this paper, we adopt the analytic 
case considered in this section to guide us in determin- 
ing which scale that is. Explicitly, we conservatively de- 
mand that that the effective beam transfer function must 
be greater than 0.15 in order to consider it to be prop- 
erly oversampled, and therefore independent of scanning 
strategy. We adopt the corresponding multipole moment 
to be ^hybrid, as defined in Section 12.31 Thus, the sym- 
metrized WMAP beam and HEALPix pixel window are 
used at scales for which the beam amplitude drops below 
0.15. 
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Fig. 4. — The ratio between the transfer functions derived in this paper and the nominal WMAP transfer functions for all DAs. Note 
that the DAs split into two main groups depending on focal plane position: The outer DAs, Kl, Kal and Ql— 2, all rise with £, whereas the 
inner DAs, Vl-2 and Wl-4, decrease with t. Note also the similarity between Wl and W4, between W2 and W3, and between VI and V2. 
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Note that if we included more years of WMAP scan- 
ning in our calculations, we could increase ^hybrid, since 
we would obtain more samples, and the WMAP point- 
ings do not perfectly coincide from year to year. This is 
computationally very expensive, and we have therefore 
chosen to limit our analysis to the range described here. 

As a direct illustration of this effect, we present in Fig- 
ure [2] the transfer function ratios for the VI WMAP DA 
with respect to the nominal WMAP transfer function, 
computed for both six months of observations and a full 
year of observations. The low and intermediate H be- 
haviour is the same for the two, but the pixel window 
effect becomes important earlier for the six month case 
than for the full year case. 

5. THE EFFECT OF ASYMMETRIC BEAMS IN WMAP 

We now present the main results obtained in this pa- 
per, namely the effective beam transfer functions for each 
WMAP DA, taking into account both the full asym- 
metric beam patterns and scanning strategy. These are 
shown in Figure [3] (red lines), and compared to the nom- 
inal WMAP transfer functions (dashed black lines). The 
vertical dotted line indicates ^hybrid for each case. 

Clearly, the differences between the two sets of results 
are relatively small, as no visual discrepancies are seen in 
this plot. However, in Figure[4]we plot the ratio between 
our transfer functions and the WMAP transfer functions 
for (. < ^hybrid, and here we do see small but significant 
differences between the two sets of results. 

First, we see that the ratios are essentially unity on the 
largest scales (smallest £'s), before they start diverging 
either towards high or low values at some characteristic 
scale. There are two exceptions to this trend, namely 
Wl and W4, which start diverging essentially already at 
very low fs. 

Next, the transfer functions split into two main groups: 
Our Kl, Kal and Ql-2 transfer functions are slightly 
higher than the corresponding WMAP functions, while 
the V and W DAs are slightly lower. 

Both of these general and qualitative remarks reflects 
the position of each DA in the WMAP focal plane (see 
Figure 6 of Jarosik et al. 2007 for an excellent visualiza- 
tion of the A side beams): Kl, Kal and Ql-2 are posi- 
tioned the furthest away from the optical axis, while VI- 
2 and Wl-4 are the closest. Similarly, Wl and W4 are 
positioned lower in elevation, and generally have more 
sub-structure, than W2 and W3. 

However, it should be emphasized that the overall dif- 
ferences are generally small, typically less than than 2% 
at £ < -^hybrid- Further, these differences are only signif- 
icant (again, with the exception of Wl and W4) in the 
intermediate- and high-^ ranges. 

To build up some intuitive understanding of the spa- 
tial variations caused by the asymmetric WMAP beams, 
we show in Figure [5] the difference between the fully 
asymmetrically convolved map and the corresponding 
map convolved with the symmetrized transfer function 
directly in harmonic space for one of the VI simulations. 
Thus, the two convolved maps have identical power spec- 
tra, but slightly different phases. The top panel shows 
the full-sky difference map with a temperature scale of 
±5/iK. The lower panels show two selected 15° x 15° re- 
gions centered on the north ecliptic pole (NEP; top row) 
and the Galactic center (GC; bottom row), respectively. 



The left column shows the actual temperature map con- 
volved with the asymmetric beam, and the right column 
shows the same differences as in the full-sky plot. 

The first striking feature seen in this map is that the 
differences are clearly larger in the ecliptic plane than 
around the ecliptic poles. This is due to the WMAP 
scanning strategy, which leads to a larger number of ob- 
servations per pixels around the poles, and also with a 
greater range of beam orientations. Next, it is difficult 
to spot any single unambiguous and well-defined corre- 
lation between the convolved and the difference maps. 
Clearly, there are similar morphological structures in the 
two, but the sign of the correlations appears to vary. 
Third, we see a clear tendency of diagonal striping in the 
GC plot, which corresponds to correlations along eclip- 
tic meridians and lines of constant latitude. (Note that 
these plots are shown in Galactic coordinates, while the 
WMAP scanning strategy is nearly azimuthally symmet- 
ric in ecliptic coordinates.) 

In the next section, we consider the impact of the 
asymmetric beams on cosmological parameters. How- 
ever, before concluding this section we make a com- 
ment concerning an outstanding issue re garding the 3- 
year W MAP power spectra first noted bv lEriksen et al.l 
(|2007t ). They pointed out the presence of a 3cr amplitude 
discrepancy between the V- and W-band power spec- 
tra (Figure 5 of Eriksen et al. 2007). Specifically, the 
V-band spectrum was biased low compared to the W- 
band spectrum between £ = 300 and 600 by ~ 80/xK^. 
'Huffe nberger et al.l ()2006[ ) later showed that - 30^K^ of 
this discrepancy could be attributed to over-estimation of 
point source power in the 3-year WMAP spectrum anal- 
ysis, and this was subsequently confirmed and corrected 
bv lHinshaw et al.l ([2007). StiU, about bO^K^ of this dif- 
ference remained, which was statistically significant at 
^ 2cr. 

lEriksen et al.l (|2007l ) proposed that this difference 
might be due to errors in the beam transfer functions 
caused by asymmetric beams. Given the new results pre- 
sented in this paper, we are now in a position to consider 
this issue more quantitatively. The relevant question 
is then whether the WMAP V-band transfer functions 
are systematically biased high compared to the W-band 
functions. At first glance, one may get this impression 
from the plots shown in Figure U) The V-band ratios 
both drop noticeably from i = 300 (decreasing nearly 
linearly from -0.2 to -0.7%), whereas W2 and and W3 are 
slightly high in the same range, at about -1-0.1 to -1-0.2%. 
On the other hand, Wl and W4 are even lower than the 
V-band functions, at -0.4 to -0.6%. The net difference is 
therefore not more than a few tenths of a percent, which 
corresponds to ^ 10/xK^ in the power spectrum. Thus, 
it is possible that this effect may contribute somewhat 
to the power spectrum discrepancy between V- and W- 
band, but it does not seem to fully explain the difference. 

6. IMPACT ON COSMOLOGICAL PARAMETERS 

We now assess the impact of asymmetric beams in 
WMAP on cosmological parameters. We do this by mod- 
ifying the co-added 5-year WMAP temperature power 
spectrum (Noha et al. 2009) pro vided w ith the WMAP 
likelihood code iDunklev et al.l (|2009l ): iKomatsu et al.l 
poos) according to the transfer function ratios shown 
in Figure IH and run CosmoMC (|Lewis fc Bridld[200l 



Convolved map Difference map 




-300 mK^b ^300 mK -5 ^K^k. ,_5 [iK 



Fig. 5. — Top panel: Difference between a VI simulation convolved with the full asymmetric beam and the same realization convolved 
with the corresponding symmetrized transfer function. The two maps have identical power spectrum but different phases. Note that 
larger differences are observed along the ecliptic plane, where the density of observations is lower than towards the ecliptic poles, and the 
cross-linking is also weaker. Bottom panels: Zoom-in on two regions, the north ecliptic pole (NEP; top row) and the Galactic center (GC; 
bottom row). Left column shows the map convolved with an asymmetric beam, and right column shows the same difference as in the top 
panel. 
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400 600 
Multipole moment, / 

Fig. 6. — Total correction to the 5-year co-added WMAP tem- 
perature power spectrum due to asymmetric beams. Note the tran- 
sition between high and low signal-to-noise weighting schemes at 
I = 500, and also the manually capped amplitude at £ > 750. The 
latter is imposed in order to be conservative in the very high-£ 
regime, where the transfer functions are sensitive to pixel window 
effects. 

to estimate the resulting parameters. Only a simple 
6-parameter ACDM model is considered in this paper. 
For comparison, we also run the code with the nominal 
WMAP spectrum as input, so that we can directly esti- 
mate the impact of asymmetric beams with everything 
else held fixed. 

Unfortunately, the individual cross-spectra for each 
pair of DAs have not yet been published by the WMAP 
team, but only the total co-added spectrum. We must 
therefore make a few approximations in order to apply 
the proper beam corrections to the full spectrum. First, 
let fj^ denote the white noise level of DA i (see Table 
[T]), li\ the transfer function estimate derived in this pa- 
per, and let /3]'^^^^ be the nominal WMAP transfer 
function. Finally, define 



if, 



0\ 



,WMAP 



1 



(13) 



to be the fractional difference between the two. 

Next, the WMAP team uses t he M ASTER pseudo- 
spectrum algorit hm (iHiv on ct al"]|2002f) for power spec- 
trum estimation ( Hinshaw et alll2003l l2007t iNolta et all 
HOO^), which quickly produces good estimates at high £'s. 
However, this method is not a maximum-likelihood esti- 
mator, and it does not yield optimal error bars. To im- 
prove on this, the WMAP applies different pixel weights 
in different multipole regions: At low £'s, where the sky 
maps are signal dominated, they apply equal weights 
to all pixels, while at high £'s, where the maps are 
noise dominated, they apply inverse noise variance pix- 
els weights. These weights are then taken into account 
when co-adding the cross-spectra obtained from all pos- 
sible DA pairs (but excluding auto-correlations). The 
transition is made at £ — 500. 

The beam-convolved (but noiseless) power spectrum 
C^'' observed by a given DA pair, i and j, may be written 
as C'^-' = Pl(3jCe, where is the true power spectrum of 
our sky, and (31 is the true transfer function for DA i. The 
noise amplitude of the same spectrum is proportional to 
^n'^nl PiPl- The inverse noise variance weight of this 



cross-spectrum is therefore approximately 



i'<j' 



(14) 



where the sum runs over all N different pairs of cross- 
spectra. (Note that this is only an approximation to 
the exact expression, because other effects also enter the 
full calculations. One important example is the sky cut, 
which couples different £ modes, and is taken into account 
through a coupling matrix. Such effects are not included 
in the analysis presented here. 

Pulling all of this together, the appropriately co-added 
power spectrum provided by WMAP should ideally read 



£ < 500 
for £ > 500 



(15) 



However, the spectrum that in fact is provided by 
WMAP is Equation [H] evaluated for Pa = /3^map^ 
which, according to our calculations, is slightly bi- 
ased. To obtain the appropriate correction factor, a, = 
Q/CWMAP^ for each £, we therefore set /Si = + 
S{) in Eguation llSl and expand to first order in Si. Doing 
this, we find that 



ai 



where u;WMAP 




{Si + S^,) for £< 500 

(Sl + Sj) for £ > 500 



,,,WMAP 



(16) 



y is the expression given in Equation [14] 
evaluated with pj^^^^ . 

This function is plotted in Figure [51 Note, however, 
that we have capped the function by hand ai £ = 750 
to be conservative, considering that our V-band trans- 
fer functions do not have support all the way to the 
maximum multipole used in the WMAP likelihood code, 

Cax = 1000. 

The results from the corresponding CosmoMC analy- 
ses are tabulated in Table [2| in terms of marginal means 
and standard deviations, and plots of the marginal dis- 
tributions are shown in Figure [Tj Here we see that there 
are small but noticeable shift in several parameters. For 
example, there is a positive shift of 0.4ct in the amplitude 
of scalar perturbations. As, and 0.3cr in the spectral in- 
dex of scalar perturbations. Although relatively modest, 
these shifts are certainly large enough that they should 
be taken properly into account. 

7. CONCLUSIONS 

This paper has two main goals. First, we wanted to 
generate a set of WMAP-like simulations that fully take 
into account the asymmetric beams and anisotropic scan- 
ning pattern of the WMAP satellite. Such simulations 
are extremely valuable for understanding the impact of 
beam asymmetries on various statistical estimators and 
models. One example of such, which indeed provided us 
with the initial motivation for studying this issue, is the 
anisotropic universe model presented bv lAckerman et al.l 
( 2007) , and late r considered in detail with resp ect to the 
WMAP data bv lGroeneboom fc Erik"senl ([2009"). The re- 
sult from that analysis was a tentative detection of vi- 
olation of rotational invariance in the early universe, or 
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TABLE 2 

COSMOLOGICAL PARAMETERS 



Parameter 


Nominal WMAP 


Corrected beams 


Shift in a 




0.0228 ± 0.0006 


0.0228 ± 0.0006 


0.1 


O fi2 


0.109 ±0.0006 


0.112 ±0.006 


0.4 


log(lOiOAs) 


3.064 ± 0.042 


3.079 ± 0.042 


0.4 


T 


0.089 ±0.017 


0.089 ±0.017 


0.00 


h 


0.722 ±0.027 


0.716 ±0.026 


-0.3 




0.965 ±0.014 


0.969 ±0.014 


0.3 



Note. — Comparison of cosmological parameters derived from 
the nominal WMAP power spectrum (second column) and from 
the power spectrum corrected for asymmetric beams (third col- 
umn). The rightmost column shows the relative shift between the 
two in units of cr. 



some other effect with similar observational signatures, 
at the 3.8(7 confidence level. It was shown that neither 
foregrounds nor correlated noise could generated this sig- 
nal, but the question of asymmetric beams was left unan- 
swered. This issue will now be revisited in an upcoming 
paper, using the simulations generated here. 

The second goal of the paper was to assess the impact 
of beam asymmetries on the WMAP power spectrum 
and cosmological parameters. We did this by comparing 
the power spectrum of the full beam convolved simula- 
tions with the power spectrum of the input realizations, 
thereby providing a direct estimate the effective beam 
transfer functions. Doing so, we found differences at the 
1 — 2% level in all differencing assemblies at intermediate 
and high £'s with respect to the nominal WMAP transfer 
functions. 

A similar analy sis was perform e d for t he 3-year WMAP 
data release by iHinshaw et al.l (|2007[ ). who approach 
the problem from an analytical point of view. How- 
ever, at that tim e only the A-side beams were available 
(jHill et al.ll2009f ). and they therefore assumed identical 
beams on both the A and B sides. With this data. 



they concluded that the impact of beam asymmetries 
was ^1% everywhere below £ — 1000 for the V- and 
W-band DAs. For comparison, we find that there is a 
~ 1% bias already a,t i — 600 for the combined co- added 
temperature power spectrum, and increasing rapidly to 
~ 1.5% a.t £ — 750. It is not unlikely that this trend 
may continue further in £, but to answer that question 
considerably more computational resources is required. 
Nevertheless, the two analyses appear to be in reason- 
able agreement with each other, especially considering 
the fact that we take into account the full beam maps of 
both the A and B sides. 

As far as cosmological parameters go, the impact of 
asymmetric beams appear to be small but noticeable. 
Specifically, we find shifts of OAa in the amplitude of 
scalar perturbations. As, and the physical density of cold 
dark matter, flcdmh^i and 0.3a in the spectral index of 
scalar perturbations, rig. While these shifts are relatively 
modest, they are of the same order of magnitude or larger 
than, say, marginalization over the Sunyaev-Zeldovich 
effect (.S pcrgcl ct al. 200^ or unresolved point sources 
(jNolta et al.ll2009l) . which indeed are taken into account. 
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One outstanding question that still remains is the im- 
pact of the anisotropic effective pixel window. As explic- 
itly demonstrated in this paper, the difference between 
the isotropized HEALPix pixel window and the actual 
WMAP VI scanning induced pixel window becomes vis- 
ible at £ ~ 900 for one year of WMAP observations. Of 
course, this is well within the noise-dominated regime for 
the WMAP data, and unlikely to have any major impact 
on cosmological results, but we believe that a proper un- 
derstanding of this issue, both with respect to WMAP 
and Planck, is warranted, and intend to revisit this issue 
in a separate study. 

The simulations described in this paper may be down- 
loaded from IKW's homepage^. 
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APPENDIX 

CONVERGENCE OF THE DIFFERENTIAL MAP MAKER 

As described in Section 12. 2i we introduce one new step to the differential map making algorithm presented by 
I Wright et al.l (|1996( l: We initialize the iterations at the exact solution of Equation [7] evaluated at low resolution, which 
in this paper is taken to be iVsido — 16, with 3072 pixels. 

To demonstrate the improvement in convergence due to this choice of initialization, we revisit the analytic case 
considered in Section [4l which compared the results from our simulation pipeline with an exact analytic case, but 
taking into account the actual WMAP scanning strategy. 

In Figure [5] we show a set of difference maps taken between the intermediate solutions produced by the differential 
map maker and the analytic and isotropic map solution. From top to bottom, the panels show the residuals after 2, 
5 and 10 iterations, and at the bottom, the final converged solutions. The left panel shows the series obtained when 
initializing the search at the low-resolution solution, while the right panel shows the series when initializing at zero. 
Convergence was achieved respectively after 67 and 123 iterations in the two cases. 

Note that the WMAP team initializes their search at the CMB dipole, which is the dominant component in their 
data set. However, this is in our setting equivalent to initializing at zero, since our simulation does not include a 
dipole. 

Taking the difference between the two final solutions, we have verified that the peak-to-peak residuals in the two 
maps are less than 0.1 /zK, of which essentially all is concentrated in a single dipole component. The solution is thus 
independent of initialization, and the only difference lies in computational speed. 

Finally, note that even though the two maps are internally indistinguishable, they are both quite different from the 
isotropic reference map. To be precise, the RMS difference between the derived maps and the isotropic reference map 
is 0.91 /LtK, with a spatial pattern similar to the overall WMAP scanning pattern. 

The cause of these residuals is once again the differences in the treatment of the effective pixel windows: The 
HEALPix pixel window is computed by uniformly averaging over the full sky, whereas the simulation pipeline takes 
into account the actual pointing directions of the satellite. Sub-pixel variations in the CMB sky therefore leads to 
significant differences in the two estimates on small scales. The effect of such pixel window variations on the 5-year 
WMAP power spectrum will be considered in a future paper. 

FAST 2D SPLINE EVALUATION 

The heart of the simulation pipeline described in Section[2]is the real-space convolution algorithm defined by Equation 
m For this operation to be computationally feasible we have to be able to evaluate the beam response function quickly 
at any position. The real beam maps, however, are provided to us in the form of two-dimensional pixelized images 
with relatively coarse resolution. It is therefore necessary to establish a fast and accurate interpolation scheme. 

We adopt a bicubic spline for this purpose, and review here one sp ecific implementation of this concept. Note that 
most of the following is standard textbook material fe.g.. ! Presil2002l) . and is included here only for easy reference. 

Suppose we are given some tabulated two-dimensional function f{x,y) over a regular grid, and want to interpolate 
at arbitrary positions (a:o,2/o) within this grid. One particularly appealing approach for doing so are by means of 
bicubic splines, which are bi-cubic polynomials in x and y, 

3 3 

P{x,y) =^^a^jx'y^ . (Bl) 

The coefficients aij are defined separately for each grid cell, and our task is to compute these given the tabulated 
function f{x, y). Note that once we have these coefficients, any spline evaluation will be very fast, since it essentially 
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Fig. 8. — Comparison of convergence of the differential map maker for two different choices of initialization. The left column shows 
the snapshots from the series obtained with initializing at a solution obtained by brute-force evaluation at low resolution, while the right 
column shows the series obtained when initializating at zero. Each plot is a difference map between the current solution for a data 
set including asymmetric beams and real scanning strategy and the corresponding map convolved with the analytic Gaussian beam and 
isotropic HEALPix pixel window. The bottom row shows the final solutions obtained in the two cases, which were obtained after 67 and 
123 iterations, respectively. These final maps are idential up to a 0.1/xK dipole. 



amounts to performing a vector-matrix-vector multiplication with a 4 x 4 matrix. 

Let us first consider a cell defined over the unit square, having corners (x,?/) — (0,0), (1,0), (0, 1) and (1, 1). (Note 
that this assumption does not imply any restriction of the problem, since any grid cell in a regular grid may be linearly 
transformed into the unit square.) Assume also that we know the function values /(x, y) and the first- and second-order 
derivatives, fx{x, y), fy{x^ y) and fxy{x, y) at all four corners. (Here subscript x denotes derivatives, fx = df /dx.) The 
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coefficients of tlie bicubic spline are tlien defined such that that both the function values and the derivatives match, 

fix,y) ^p{x,y) (B2) 
fx{x,y) =Px{x,y) (B3) 
fy{x,y) ^Py{x,y) (B4) 
fxy{x,y) = Pxy{x,y), (B5) 

at all four corners. With four equations for each of four corners, there is a total of 16 independent equations for the 
16 independent spline coefficients, a^ , and the spline is therefore uniquely defined. 
Writing out Equation IB5I explicitly, we obtain the following set of linear equations. 



f(0,0) = 


n(0,0) = 


/(1,0) = 


p(i,o) = 


/(o,i) = 


P(0,1) = 


/(!,!) = 


= 


/x(0,0) = 


Px(0,0) = 


/.(1,0) = 


P.(1,0) - 


/x(0,l) = 


Px(0,l) = 






fy{0,0) = 


p,(0,0) = 


/.(1,0) = 


Py(l,0) = 


/.(0,1) = 


P?y(0,l) = 




Py(l,l) = 


^,(0,0) : 


= Pxy(0,0) 


/.y(l,0) = 




fxy{0,l) 


= Pa;y(0, 1) 


fxy{l,l) 
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The remaining problem is then to estimate the first- and second-order derivatives at all grid points, and several 
approaches may be used for this. We adopt a spline based method for this step as well. 

First, we compute a standard one-dimensional natural spline along all x and y coordinate lines, using standard 
methods. The result from this operation is the set of all pure second-order derivatives, fxx{x, y) and fyy{x, y), at each 
grid point. 

Because a one-dimensional spline is also a simple cubic polynomial, and therefore only has four free coefficients, 
it is sufficient to know the function values and second-order derivatives at all grid corners to uniquely specify every 
coefficien t. As a cons equence of this, the first-order derivative along the x-direction at grid point {xi,yj) is uniquely 
given by (|Pressll2002n 

fx{x.,yj) = fi^^+^^y^)^- f(^^^y^) _ y^.)- (BIO) 

Here hx = x^+i — Xi is the x-direction grid cell size for the current cell. An equivalent expression obviously holds for 
the y-direction derivatives. 

Finally, to estimate the second-order cross-derivatives, fxy{x, y), the above process is repeated such that ^/-derivatives 
are computed from fx{x,y) splines for all y-direction coordinate lines. Thus, all required derivatives may be obtained 
by performing m + 2n one-dimensional spline computations, where m is the number of grid cells in x-direction, and n 
is the number of grid cells in y-direction. 
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